function [pM,pS,pS_2w,pi,V,Vtp1] = solve_problem(states,fixedparms,xadata,spec,beliefparms) 

% solve dynamic problem recursively

    % set seed so simulate same eps on each iteration
    aa=0;
    rng(1)
    
    % set an indicator for which states the daughter is still in school
    inschool=(states(:,3)==1);
    
    % set up terminal utility matrix
    %-------------------------------
    % set up the utility lookup matrix. Gives the probability of choosing
    % the terminal action in every decision state when faced with the low
    % quality offer (column 1) and the high quality offer (column 2)
    u=utility_f(states,xadata,fixedparms,spec);   

    % set up pi matrix 
    %------------------
    pi = makepi(states,fixedparms,beliefparms);

    % pi override (must be commented out for everything but running
    % counterfactuals
    %------------
%     global pi_counter
%     
%     if fixedparms.ind_goodgirl==1 
%         pi=pi_counter(:,1:2)
%     end
%     
%     if fixedparms.ind_goodgirl==0 
%        pi=pi_counter(:,3:4)
%     end
%     pi

    % set up choice probability matrix - defined for each state
    %----------------------------------
    pM=NaN(size(states,1),2);
    pS=NaN(size(states,1),2);
    % choosing school in the two way choices if don't get marriage offer
    pS_2w=NaN(size(states,1),1);
    
    % vector for value functions - defined as the value of being unmarried at the
    % beginning of period t before know type of marriage offer that will
    % get
    %----------------------------
    V=NaN(size(states,1),1);
    
    % V(t+1) -- columns: one conditional on choosing home, one on school
    Vtp1=NaN(size(states,1),2);
    
    % fill in last period
    Vtp1(states(:,1)==22,1)=beliefparms.V0.*u(states(:,1)==22,1)+beliefparms.VTp1;

    
    % begin with recursive problem for all two-way choices in states when are out of school
    %-------------------------------
    for a=22:-1:13
        
            pM(states(:,1)==a & inschool==0,1)=normcdf((u(states(:,1)==a & inschool==0,1)- Vtp1(states(:,1)==a & inschool==0,1))./(fixedparms.disc.^(a-12).*(beliefparms.sigma_ea.^2+beliefparms.sigma_ea.^2).^0.5));
            pM(states(:,1)==a & inschool==0,2)=normcdf((u(states(:,1)==a & inschool==0,2)- Vtp1(states(:,1)==a & inschool==0,1))./(fixedparms.disc.^(a-12).*(beliefparms.sigma_ea.^2+beliefparms.sigma_ea.^2).^0.5));
        
            % use these calculated probabilities to calculate value function for
            % period t
            %-------------------------------------------------------------------
            % first, calculate conditional expectations for the epsilons
            z0=(Vtp1(states(:,1)==a & inschool==0,1)-u(states(:,1)==a & inschool==0,1))./(fixedparms.disc^(a-12));
            z1=(Vtp1(states(:,1)==a & inschool==0,1)-u(states(:,1)==a & inschool==0,2))./(fixedparms.disc^(a-12));
        
            E_eps_a_0=(fixedparms.disc^(a-12)).*(1./(2.^0.5)).*beliefparms.sigma_ea.*normpdf(z0./(2.*beliefparms.sigma_ea.^2).^0.5)./(1-normcdf((z0./(2.*beliefparms.sigma_ea.^2).^0.5)));
            E_eps_a_1=(fixedparms.disc^(a-12)).*(1./(2.^0.5)).*beliefparms.sigma_ea.*normpdf(z1./(2.*beliefparms.sigma_ea.^2).^0.5)./(1-normcdf((z1./(2.*beliefparms.sigma_ea.^2).^0.5)));
            E_eps_r_0=(fixedparms.disc^(a-12)).*(1./(2.^0.5)).*beliefparms.sigma_ea.*normpdf(-z0./(2.*beliefparms.sigma_ea.^2).^0.5)./(1-normcdf((-z0./(2.*beliefparms.sigma_ea.^2).^0.5)));
            E_eps_r_1=(fixedparms.disc^(a-12)).*(1./(2.^0.5)).*beliefparms.sigma_ea.*normpdf(-z1./(2.*beliefparms.sigma_ea.^2).^0.5)./(1-normcdf((-z1./(2.*beliefparms.sigma_ea.^2).^0.5)));
             
            % replace with max if NaN
            E_eps_a_0(isnan(E_eps_a_0))=2.*beliefparms.sigma_ea;
            E_eps_a_1(isnan(E_eps_a_1))=2.*beliefparms.sigma_ea;
            E_eps_r_0(isnan(E_eps_r_0))=2.*beliefparms.sigma_ea;
            E_eps_r_1(isnan(E_eps_r_1))=2.*beliefparms.sigma_ea;
        
            E_eps_a_0(E_eps_a_0>2.*beliefparms.sigma_ea)=2.*beliefparms.sigma_ea;
            E_eps_a_1(E_eps_a_1>2.*beliefparms.sigma_ea)=2.*beliefparms.sigma_ea;
            E_eps_r_0(E_eps_r_0>2.*beliefparms.sigma_ea)=2.*beliefparms.sigma_ea;
            E_eps_r_1(E_eps_r_1>2.*beliefparms.sigma_ea)=2.*beliefparms.sigma_ea;
        
            pi0=pi(states(:,1)==a & inschool==0,1);
            pi1=pi(states(:,1)==a & inschool==0,2);
            
            u0=u(states(:,1)==a & inschool==0,1);
            u1=u(states(:,1)==a & inschool==0,2);
            
           % acceptance probabilitie
            phat0=pM(states(:,1)==a & inschool==0,1);
            phat1=pM(states(:,1)==a & inschool==0,2);
        
            %V(t)
            V(states(:,1)==a & inschool==0)=pi0.*phat0.*(u0+E_eps_a_0)...
                + pi1.*phat1.*(u1+E_eps_a_1)...
                + pi0.*(1-phat0).*(Vtp1(states(:,1)==a & inschool==0,1)+E_eps_r_0)...
                + pi1.*(1-phat1).*(Vtp1(states(:,1)==a & inschool==0,1)+E_eps_r_1)...
                +(1-pi0-pi1).*Vtp1(states(:,1)==a & inschool==0,1);         
   
        
            % fill in Vtp1 for a-1
            %----------------------
            for e=7:13
                if size(Vtp1(states(:,1)==a-1 & states(:,2)==e &  inschool==0,1),1)>0
                    Vtp1(states(:,1)==a-1 & states(:,2)==e &  inschool==0,1)=V(states(:,1)==a & states(:,2)==e &  inschool==0,:);
                    end
                end
            
    end
    


     
    % now go onto the three way choices, beginning at age 18 
    %--------------------------------------------------------
    % at age 18, when still in school, today's choice is three-way but tomorrow's choice is
    % twoway. VTp1(18, inschool| Ed)=V(20,13)
    Vtp1(states(:,1)==18 & states(:,4)==13,2)=V(states(:,1)==19 & states(:,2)==13);
    Vtp1(states(:,1)==18 & states(:,4)==13,1)=V(states(:,1)==19 & states(:,2)==12);
   
    
    % loop recursively through all the three way choices
    %---------------------------------------------------

   
    % simulate trivariate distribution of preference shocks - we will
        % need these to estimate the expected value of the error
        % conditional on the action being optimal
    
        rng(1)
        mu=[0 0 0];
        cov=[beliefparms.sigma_ea^2 0 0; 0 beliefparms.sigma_ea^2 0; 0 0 beliefparms.sigma_ea^2];
        eps=mvnrnd(mu,cov,10000);
        eM=eps(:,1);
        eH=eps(:,2);
        eS=eps(:,3);

    for a=18:-1:13
        
        % create indicator to mark which choice we are interested in
        ind=(states(:,1)==a & inschool==1);
    
        % given Vtp1, estimate choice probabilities 
        % this is the covariance of the difference between pref shocks
        cov=[2*beliefparms.sigma_ea^2 beliefparms.sigma_ea^2; beliefparms.sigma_ea^2 2*beliefparms.sigma_ea^2];
        mu=[0 0];
         
        pM(ind==1,1)=mvncdf([u(ind==1,1)-Vtp1(ind,1) u(ind==1,1)-Vtp1(ind,2)], mu , (fixedparms.disc.^(a-12))^2.*cov);
        pM(ind==1,2)=mvncdf([u(ind==1,2)-Vtp1(ind,1) u(ind==1,2)-Vtp1(ind,2)], mu , (fixedparms.disc.^(a-12))^2.*cov);
        
        pS(ind==1,1)=mvncdf([Vtp1(ind,2)-Vtp1(ind,1) Vtp1(ind,2)-u(ind==1,1)], mu , (fixedparms.disc.^(a-12))^2.*cov);
        pS(ind==1,2)=mvncdf([Vtp1(ind,2)-Vtp1(ind,1) Vtp1(ind,2)-u(ind==1,2)], mu , (fixedparms.disc.^(a-12))^2.*cov);
    
        % two way choice probs if don't get an offer
        pS_2w(ind==1)=normcdf((Vtp1(ind,2)-Vtp1(ind,1))./((fixedparms.disc.^(a-12)).*beliefparms.sigma_ea.*2^0.5));
              
           
    
        % using these probabilities, form the value function for three way
        % choices
        %------------------------
        % utilities associated with accepting marriage offer
        u0=u(ind==1,1);
        u1=u(ind==1,2);
        
        % probability of getting each type of marriage offer
        pi0=pi(ind==1,1);
        pi1=pi(ind==1,2);
        
        % choice probabilities conditional on type of marriage offer
        phatM0=pM(ind==1,1);
        phatM1=pM(ind==1,2);
        phatS0=pS(ind==1,1);
        phatS1=pS(ind==1,2);
        phatS_2w=pS_2w(ind==1);
        
        % next period's value function conditional on choice
        Vtp1S=Vtp1(ind==1,2);
        Vtp1H=Vtp1(ind==1,1);
        
        % use simulated variables to solve for E(eps) in 3w choices
        E_eps_3wM_0=(fixedparms.disc^(a-12)).*mean(eM(eH-eM<(u0-Vtp1H)./(fixedparms.disc^(a-12)) & eS-eM<(u0-Vtp1S)./(fixedparms.disc^(a-12))));
        E_eps_3wH_0=(fixedparms.disc^(a-12)).*mean(eH(eM-eH<(Vtp1H-u0)./(fixedparms.disc^(a-12)) & eS-eH<(Vtp1H-Vtp1S)./(fixedparms.disc^(a-12))));
        E_eps_3wS_0=(fixedparms.disc^(a-12)).*mean(eS(eM-eS<(Vtp1S-u0)./(fixedparms.disc^(a-12)) & eH-eS<(Vtp1S-Vtp1H)./(fixedparms.disc^(a-12))));
        
        E_eps_3wM_1=(fixedparms.disc^(a-12)).*mean(eM(eH-eM<(u1-Vtp1H)./(fixedparms.disc^(a-12)) & eS-eM<(u1-Vtp1S)./(fixedparms.disc^(a-12))));
        E_eps_3wH_1=(fixedparms.disc^(a-12)).*mean(eH(eM-eH<(Vtp1H-u1)./(fixedparms.disc^(a-12)) & eS-eH<(Vtp1H-Vtp1S)./(fixedparms.disc^(a-12))));
        E_eps_3wS_1=(fixedparms.disc^(a-12)).*mean(eS(eM-eS<(Vtp1S-u1)./(fixedparms.disc^(a-12)) & eH-eS<(Vtp1S-Vtp1H)./(fixedparms.disc^(a-12))));
        
        E_eps_2wS=(fixedparms.disc^(a-12)).*mean(eS(eH-eS<((Vtp1S-Vtp1H)./(fixedparms.disc^(a-12)))));
        E_eps_2wH=(fixedparms.disc^(a-12)).*mean(eS(eS-eH<((Vtp1H-Vtp1S)./(fixedparms.disc^(a-12)))));
    
        % cap the expected value of the preference shocks at 2*sigma
        % (discounted appropriately) 

        E_eps_3wM_0(isnan(E_eps_3wM_0))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_3wH_0(isnan(E_eps_3wH_0))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_3wS_0(isnan(E_eps_3wS_0))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_3wM_1(isnan(E_eps_3wM_1))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_3wH_1(isnan(E_eps_3wH_1))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_3wS_1(isnan(E_eps_3wS_1))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_2wS(isnan(E_eps_2wS))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_2wH(isnan(E_eps_2wH))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
    
        E_eps_3wM_0((fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea<(E_eps_3wM_0))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_3wH_0((fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea<(E_eps_3wH_0))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_3wS_0((fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea<(E_eps_3wS_0))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_3wM_1((fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea<(E_eps_3wM_1))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_3wH_1((fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea<(E_eps_3wH_1))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_3wS_1((fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea<(E_eps_3wS_1))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_2wS((fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea<(E_eps_2wS))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 
        E_eps_2wH((fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea<(E_eps_2wH))=(fixedparms.disc^(a-12)).*2.*beliefparms.sigma_ea; 

    
        % solve value function  
        V(ind==1)=pi0.*(phatM0.*(u0+E_eps_3wM_0)+phatS0.*(Vtp1S+E_eps_3wS_0)+(1-phatM0-phatS0).*(Vtp1H+E_eps_3wH_0))... 
            + pi1.*(phatM1.*(u1+E_eps_3wM_1)+phatS1.*(Vtp1S+E_eps_3wS_1)+(1-phatM1-phatS1).*(Vtp1H+E_eps_3wH_1))... 
            +(1-pi0-pi1).*(phatS_2w.*(Vtp1S+E_eps_2wS)+(1-phatS_2w).*(Vtp1H+E_eps_2wH));
        
        % role this value function back to t-1
        if a>13
            % if choose to stay in school
            Vtp1(states(:,1)==a-1 & states(:,3)==1 ,2)=V(ind==1);

            % if choose to drop out of school
            edcurrent=states(states(:,1)==a-1 & states(:,3)==1 ,2);
            Vtp1(states(:,1)==a-1 & states(:,3)==1 ,1)=V(states(:,1)==a & states(:,2)==edcurrent);
        end

    end

end